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ABSTRACT 

NGC 4013 is a distinctly warped galaxy with evidence of disk-halo activity. Through deep H I 
observations and modeling we confirm that the HI disk is thin (central exponential scale height of 
with an upper limit of 4” or 280 pc), but flaring. We detect a vertical gradient in rotation velocity 
(lag), which shallows radially from a value of —35^28 km s -1 kpc -1 at 1.4’ (5.8 kpc), to a value of zero 
near R25 (11.2 kpc). Over much of this radial range, the lag is relatively steep. Both the steepness 
and the radial shallowing are consistent with recent determinations for a number of edge-ons, which 
have been difficult to explain. We briefly consider the lag measured in NGC 4013 in the context of 
this larger sample and theoretical models, further illuminating disk-halo flows. 


1. INTRODUCTION 

To fully comprehend star formation (SF) processes, 
we must understand the properties and origins of the 
gas that fuels them. If galaxies are to avoid exhausting 
their star forming fuel in a short time compared to their 
age, then there must exist some means of gaining mate¬ 
rial, possi bly through _stellar mass loss as suggested by 
iLeitner fc Kravtsov! (1201 11 1. or through the accretion of 
primordial neutral and ionized gas. Thus, extra-planar 
layers are of great interest as they are the interface be¬ 
tween the main disk and the intergalactic medium (IGM) 
and all material accreted onto the midplane to fuel SF 
must first pass through them. 

With the exception of neutral hydrogen Hi, a 
connection has been made between the presence of 
extra-planar compo nents and SF both on local and 
global scales (e.g. iPutman et al.l 120121 and references 
therein). In particular, connec tions h ave been made 
between SF and hot gas (e.g . jTiillmann etahl 120001 
iHodges-Kluck fc Bregmanl I2013L iLi fc Wand 1201 3l ) . rel¬ 

ativistic particles (e. g. Ilrwin et al.1 I1999D dust (e.g. 
iHowk fe Savagel I1999D. and diffuse ioni zed gas (DIG) 
(e.g. lRandlll996tlRossa fe Dettmad 120031) . Extra-planar 
H 1 has been shown to have some c onnection to disk- 
hale^ flo ws wi thin some galaxies (e.g. Ilrwi 3 (j!994T > and 
Boomsm a et all (1200511. but e vidence for Hi accretion 
[e.g. HVCs in the Milky Way dWakker fe van Woerdenl 
199711 . counterrotatin g clouds in NGC 891 seen by 
Oosterloo et al.ll2007j| exists, leaving the relative contri¬ 
butions of each process ambiguous until further study. 

A powerful method to determine the origins of extra- 
planar layers is to study their kinematics. In particular, 
decreases in rotation speed with height above the mid- 
plan e (lags) have come to the forefront in r ecent years 
(e.g. lOosterloo et aLll2007l iHeald et al.ll2007ll leading to 
various models attempting to explain them in terms of 
disk-halo flows or an interaction between disk-halo cy¬ 
cling gas and accreting gas. Purely ballistic models of 
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disk halo flow (e.g. iCollins et all I2002L iBarnabe et al.l 
120061 IFraternali fc Binnevl 120061) consider only gravita¬ 
tional effects: Material is ejected from the midplane in 
disk-halo flows and then moves outward radially due to 
its being farther from the gravitational potential. To 
conserve angular momentum, the rotational velocity of 
the material decreases, creating a lag with height above 
the midplane. A clear shortcoming of such models is that 
they greatly under- predict lag m agnitudes by up to an 
order of magnitude (IHeald et al.ll2007 ). Thus, additional 
processes may be at work such as pressure gradients or 
magnetic tension (Benjamin 2002, 2012). However, sim¬ 
ulations involving these mechanisms have yet to be de¬ 
veloped in earnest, with future progress made possible 
through observational constraints, in particular those in¬ 
volving lags. 

Pressure gradients and magnetic tension aside, some 
theoretical simulations have gone beyond purely ballis¬ 
tic effects and disk -halo flows cons i dered in isolation. 
Those presented by iMarinacci et al.l (1201 If) involve mo¬ 
mentum and heat exchanges occurring in disk-halo flows. 
When cool clouds are launched from the disk into a 
hot, initially static corona built up by accretion, gas is 
stripped from the clouds, which streams in the cloud’s 
wake, mixing with the hot coronal gas. The mixed gas 
forms additional HI clouds before falling to the mid¬ 
plane. This mixing results in substantially more efficient 
cooling of the hot coronal gas than_is descri bed in the 
earlier, related model of IFraternali fc Binnev (200 8). in 
which hot coronal gas experienced cooling times ~100 
times greater than the expected flow time of galactic 
fountain clouds, thus r esulting in most of the hot gas 
remaining in the halo. IMarinacci et all (120111 1 demon¬ 
strate that this process would contribute approximately 
—7 km s -1 kpc -1 to the deceleration of clouds, amount¬ 
ing to 50% of the total —15 km s -1 kpc -1 observed 
in NGC 891. If one also considers ballistic effects, this 
momentum exchange would nearly eliminate discrepan¬ 
cies between purely ballistic models and observed lags 
in NGC 891 and the Milky Way. However, lags have 
been shown to di f fer substantially among galaxies (e.g. 
Heald et al.1120071 iSancisi et a.l.ll20(M iZschaechner et all 
20121 iKamphuis et al.l 120131) . Thus, additional galaxies 
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must be observed and modeled to comparable sensitiv¬ 
ity and detail. The Westerbork Hydrogen Accret ion in 
LOcal GAlaxieS (HALOGAS) survey (|Heald et a.l.ll 20 1 ll ) 
has sought to remedy this situation through extensive HI 
observations of a carefully selected sample of nearby spi¬ 
ral galaxies. 

Not only are lag magnitudes of interest, but variations 
in lag magnitude with radius are expected based on con¬ 
servation of angular momentum. From ballistic consider¬ 
ations only, in a disk-halo flow, material at smaller radii 
experiences larger relative changes in the gravitational 
potential when ejected from the midplane to a given 
height (and thus a gr e ater lag) than material at large 
radii. lOosterloo et al.1 (120(171) . Zsc haechner et al. (2011, 
20121. and iKamphuis et al.l (|2011l) measure a shallowing 
of the lags in several nearby galaxies, w ith a current sum¬ 
mary given in lZschaechner et al.l ( 2015 1. A general trend 
noted in that paper is that the lags begin to shallow near 
a radius of O. 5 R 25 , and reach their shallowest point near 
R 25 itself. Considering the mismatch between observed 
lag magnitudes and ballistic models, it is likely that the 
origin of the shallowing is more complicated as well. 

This work builds upon the HALOGAS sample and 
models, as well a s supp lemental galaxies presented in 
I Zschaechner et ahl (120151) . To this existing collection of 
modeled galaxies, we add NGC 4013, a nearby, well- 
studied, edge-on spiral galaxy with multiple extended 
extra-planar components. 

1 . 1 . NGC 4013 

NGC 4013 is a member of the NGC 4151 group 
containin g 16 galaxies, but is itself rather isolated 
dGiuricin et al.l I2000D . Arguably, its most noteworthy 
f eature is its s ubstan tial warp originally seen in H 1 by 
IBottema et al.l (119871) . Here we present detailed H 1 ob¬ 
servations and tilted-ring models of this galaxy. A brief 
summary of NGC 4013’s observational properties is pro¬ 
vided in Table |T| Before presenting our work, we will 
review some of the key findings throughout the litera- 
tu re con cerning NGC 4013. 

iRandl ( 19964 presents an analysis of the DIG in nine 
nearby, edge -011 spiral galaxies, including NGC 4013. In 
that work, four extra-planar filaments are found to ex¬ 
tend 2.5 kpc from the disk, defining a nearly “H” or 
“X” shape. These filaments could be related to out¬ 
flows fueled by SN activity near the center. The extra- 
planar DIG is not nearly as prominent as in NGC 891 
or NGC 5775, which have higher SFRs. The pres- 
ence of extended extr a-planar DIG is confirmed by 
iMiller fe Veilleuxl (I2003D , who also see indications of gra- 
dients in rotation speed with height above the midplane. 
IR.ueff et ahl (2013 1 observe the DIG, but also find extra- 
planar dust up to approximately 2 kpc above the disk. 
The dust is highly structured and filamentary - far more 
so than the DIG. They conclude that direct evidence for 
a connection between the extra-planar DIG and extra- 
planar dust is lacking, suggesting the gas associated with 
t he dust is either at o mic o r molecular. 

IVerstappen et al.l ( 2013 1 observe the dust in seven 
nearby, edge-on spiral galaxies using Herschel. Accord¬ 
ing to their analysis, NGC 4013 is the only galaxy within 
their sample to display substantial extra-planar dust. 

All of these results suggest disk-halo cycling at a level 
above those of galaxies such as NGC 4244 and NGC 4565, 


but low er than that of NGC 891. 

iGarcia-Burillo et all (|l999l) present CO observations 
of NGC 4013, emphasizing its box-shaped bulge that is 
likely due to a bar. They note extra-planar extensions 
showing so me c on nect ion to the extra-planar DIG de¬ 
scribed bv IRandl (119961) . Additionally, they find a 130 
km s -1 spread in velocities near the center, which they 
n ote has no corresp onding H 1 feature. 

iMartinez-Delgado et al.l ( 2009 1 discovered a tidal stel¬ 
lar stream extending 3 kpc above the plane of the disk 
in NGC 4013, possibly indicating a minor merger has 
taken place approximately 3 Gyr ago. Such a scenario 
was previously unexpected, as NGC 4013 appears to be 
relatively isolated and undisturbed. 

IR.ueff et al.l ( 2014 ) observe an extra-planar H 11 region 
in NGC 4013 that is 850 pc above the midplane and 
has a substantially lower metallicity than the main disk. 
They interpret this as a mixing of material in the gas 
that formed the associated stars. 

Perhaps most relevant to the work presented h ere, H 1 
observ ations and modeling were presented in Bottemal 
(I1995D . Those observations were done using the WSRT, 
achieving a spatial resolution of 19” x 13”, velocity reso¬ 
lution of 33.3 km s - , and a single channel rms noise of 
0.3 nrjy beam -1 . As will be shown in the next section, 
ou r data imp rove upon these numbers. 

(Bottgnyi(l996|) presented tilted-ring models of NGC 
4013 using the aforementioned WSRT data. Those mod¬ 
els include a central inclination of 90°, and a central po¬ 
sition angle of 66 °. A substantial warp is also modeled. 
The rotation curve peaks at 195 km s -1 and decreases 
to 165 km s -1 at large radii. Naturally, the models pre¬ 
sented in this paper will also constrain these parame¬ 
ters, which we derive independently but compare in § [7] 
Additionally, the models presented here focus largely on 
t he extra- pl ana r dynamics, which were not considered in 
( Bottema 199 6;). 

2. OBSERVATIONS & DATA REDUCTION 

Observations of NGC 4013 were completed with the 
VLA in B and C configurations. Data reduction was 
performed in AIPS using standard spectral line meth¬ 
ods. Self-calibration was performed after combining all 
tracks, but before continuum subtraction, and led to sub¬ 
stantial improvement. Although the velocity resolution 
of these data is 6.7 km s -1 , due to incomplete (not all of 
the awarded time was observed) observations, we aver¬ 
age three channels to obtain our desired rms noise level 
for the final cube. A variety of weighting schemes were 
used to make cubes, with the final, full resolution cleaned 
cube using Briggs weighting and a robust parameter of 
3 after attempting a range of robust parameters to find 
the optimal balance between resolution and sensitivity. 
A low spatial resolution cube (used only for the lowest 
contour in Figure [T| uses a robust parameter of 2 with 
an outer uv taper of 7 kilo-lambda. Observation dates 
and cube parameters are given in Table [2] The primary 
beam correction was performed in AIPS. 

3. THE DATA 

Immediately clear from the zerotli-moment map (Fig¬ 
ure [ 1 ]) is the prominent, nearly symmetric warp in NGC 
4013. Hints of a line of sight warp component, flare, or 
possibly a thickened H 1 layer may also be seen extend- 
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TABLE 1 

NGC 4013 Parameters 


Parameter 

Value 

Reference 

Distance (Mpc) 

14.6 a 

Nasa Extragalactic Database 

Systemic velocity (km s —1 ) 

835 

Bottema (1995) 

Inclination 

90° 

This work 

Morphological Type 

SBc 

^ull^^^ah (4988) 

SFR (Mgyr- 1 ) 

0.48 

Weigert et al. (2015, submitted to A J) 

Kinematic Center a (J2000.0) 

llh 58m 31.6s 

This work 

Kinematic Center 5 (J2000.0) 

43d 56m 52.2s 

This work 

D 25 (kpc) 

22.3 

de Vaucouleurs et al. (1991) 

Total Atomic Gas Mass (1O 9 M0) 

2.1 c 

This work 


a Distance is the median value of distances found on the NED database, excluding those obtained 
using the Tully-Fisher relation to be consistent with HALOGAS. 

k Includes neutral He via a multiplying factor of 1.36. Value obtained using single dish data at our 
adopted distance is 2.4x10 9 Mq dSpringob et aL||2005f) . WSRT value obtained by IBottemal d 1995T) 
using our adopted distance is 2.3x10 9 Mq. 


TABLE 2 

Observational and Instrumental Parameters 


Parameter 

Value 

Observation Dates — C Config. 

2010 Nov 16 
2010 Nov 20 
2010 Nov 28 

B Configuration 

2011 Mar 21-23 
2011 Apr 02 
2011 Apr 22 
2011 May 02 
2011 May 04 
2011 May 09 

Total Time on source - C Config. (hours) 

9 a 

Total Time on source - B Config. (hours) 

9.5 b 

Pointing Center 

llh 58m 31.30s 
43d 56m 48.00s 

Number of channels 

256 

Velocity Resolution 

20.1 km s -lc 

Primary Cube Beam Size 

12.7x9.91” 
900x700 pc 

PA 

4.92° 

RMS Noise — la 1 chan. 

0.18 mjy beam -1 

Secondary Cube Beam Size 

23.4x19.7” 
1700x1400 pc 

PA 

84.37° 

RMS Noise — la 1 chan. 

0.23 mjy beam -1 


a Although nine total hours were observed, approximately three of 
these hours had substantial RFI and were flagged heavily, 
k Approximately two hours of these data were omitted due to poor 
quality, resulting in 7.5 hours used in the final cubes. 
c Averaged from 6.7 km s —1 to 20.1 km s -1 for modeling in order to 
improve SNR. 


ing to an apparent height of ~0.7’ (3 kpc). Thickening 
of the disk largely to the southwest strongly suggests 
evidence for a line of sight warp component, although 
the other possibilities (or a combination of them) can¬ 
not yet be ruled out. It is unlikely that the warp is 
oriented across the line of sight - its more probable that 
the warp is indeed comprised of both across and along 
line of sight components. Further investigation through 
modeling (§ a is required to determine which of these 
components are actually present in NGC 4013. The final 
constraints on these components will be given in § [TJ 
Channel maps (Figure [2]) again show the prominent 
warp in NGC 4013. Evidence for a line of sight warp 
component may be seen by the splitting in channels cor- 
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FlG. 1.— A zeroth-moment map of NGC 4013. The full resolu¬ 
tion cube is represented in black, while an outline of the smoothed 
cube is in gray. Contour levels for the full resolution cube begin at 
2.6x10 20 cm -2 and increase by factors of two. The gray contour 
is at 9.5x10 19 cm -2 . Black lines represent the slice locations of bv 
diagrams in Figure m The beams are shown in the lower left-hand 
corner (white is full resolution, black is the smoothed cube). Note the 
prominent warp. Also note the apparent thickness of the warped re¬ 
gions, possibly due to a flare, thick disk, or a warp component along 
the line of sight. 


responding to 736 km s _1 and 756 km s — L It is clear 
that the disk is very thin near the center. Note that more 
than half of the total radial extent is clearly affected by 
the warp. The component of the warp across the line 
of sight extends a substantial 3’ (12.8 kpc) above the 
midplane at its highest point. 

The total mass we derive from our observations (Ta¬ 
ble [TJ is cl ose to that derived from single-dish data 
dSpringob et al.1120051 1 , which implies that there is no 
substantial flux missing. A first assessment going by ap¬ 
parent extent only indicates that 32% of the total HI 
mass is at a height greater than 1 kpc (14”). If only 
material within a radius of 7.9 kpc (2.2’), thus avoiding 
the clearly warped regions, this percentage amounts to 
18%. Again, we emphasize that this is only the HI with 
an apparent extent above 1 kpc. Naturally if a line of 
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sight warp component is present, much of this gas will 
be closer to the midplane. 

We now consider tilted-ring models created to fit these 
data. 

4. THE MODELS 

For tilted-ring modeling, in which we create model 
galaxies comprised of a series of concentric rings for 
which parameters such as the rotational velocity, incli¬ 
nation, position angle, etc. may be s pecified , we use the 
Tilted Ring Fitting Code ([TiRiFiG. I.Tozsa et, al.l 12007 ). 
Through tilted-ring modeling, features such as warps 
(via changes in inclination and position angle) and flares 
(vial radial increases in scale height) may be constrained. 
In addition to capabilities included in previous tilted- 
ring modeling software, TiRiFiC allows for the fitting of 
asymmetries and localized features, as well as lags, which 
we utilize in our models of NGC 4013. 

The approaching and receding halves of NGC 4013 
are modeled separately as it is clear upon initial inspec¬ 
tion that their morphologies differ greatly, and fitting 
them together could lead to unnecessary complexity in 
the models. Additional asymmetries almost certainly ex¬ 
ist, but due to projection effects, these are not obvious 
upon inspection. Thus, we only allow for asymmetry be¬ 
tween the two halves. The initial estimate for the incli¬ 
nation was 90° and the position angle (PA) for the non- 
warped disk is 65°, both determined by eye (with only 
very slight modifications to the central PA during later 
stages of the modeling and no modification to the cen¬ 
tral inclination). The initial rotation curve was estimated 
based on matching the velocities on the terminal sides of 
the data and model lv (major-axis position-velocity) di¬ 
agrams. Subsequent modifications to the rotation curve 
were made in the same manner until optimal fits were 
reached. The initial surface brightness distribution was 
set by examining the flux distribution along the midplane 
by eye. The surface brightness distribution was altered 
only slightly during the modeling, with the final distribu¬ 
tion optimized for our final preferred model, and imposed 
on all other models. The reasoning behind imposing the 
surface brightness distribution from the final model onto 
the other models is 1) any differences between models 
were minute and could be attributed to degeneracies be¬ 
tween the rotation curve and surface brightness, and 2) 
to maintain a clear and consistent picture of the changes 
in the models due to each additional component. All 
models presented here share the same surface brightness 
distribution, velocity dispersion (15 km s -1 ), systemic 
velocity (835 km s -1 ), central position (llh 58m 31.6s, 
43° 56’ 52.2”), central inclination (90°), radial position 
angle dependence, and rotation curve. The only excep¬ 
tion is that the models including a lag have a rotation 
curve that is 10 km s^ 1 higher in the approaching half 
than those that do not. The parameters not given here 
are shown in Figure [3] for final models. 

4.1. Individual Models 

We now address individual models and their defining 
characteristics, starting with a simple one-component 
model with a warp component across the line of sight 
(1C in figures). 

As stated previously, it is clear that there exists a warp 
component in NGC 4013 that lies across the line of sight. 


Thus, we model it with a single disk having an exponen¬ 
tial scale height optimally found to be 7” ±0.5 (500±35 
pc). To this model we add a PA warp for which param¬ 
eters are given in Figure [3] This model fails to match 
the vertical profile [summed over ±60” (4.3 kpc) from 
the center in order to avoid the warped outer regions as 
much as possible (Figure [3])]. The width of the verti¬ 
cal profile is overestimated, yet the model cannot repro¬ 
duce the spread near the systemic side of bv (minor-axis 
position-velocity) diagrams (Figure [5]). Thus, the model 
is too wide on the terminal side, but too narrow on the 
systemic. Similarly, as may be seen in channel maps 
shown in Figure [6l the disk is too thick at small radii, 
but too narrow at large radii. The splitting of most of 
the data bv diagrams near the systemic side is indicative 
of a warp component along the line of sight, which we 
now address. 

We add a line of sight warp component (“W” model 
in figures, inclination parameters given in Figure [3]) that 
begins near the radial onset of the other warp component 
(~1.2’ or 10 kpc), but is roughly a quarter the magni¬ 
tude. The relative contributions of each component to 
the total warp indicate that, in reality, the warp is ori¬ 
ented at an angle approximately 70 ° from the line of 
sight (of course, there is a degeneracy due to projection 
effects in terms of the near and far sides). The scale 
height of the layer is reduced to 5” (350 pc) to accom¬ 
modate for the effects of the warp in this model. One can 
see clear improvements to the bv diagrams and channel 
maps, but can also see that the emission near the sys¬ 
temic side is not entirely matched in the former, nor are 
the outer edges of the latter. The warp also improves 
the fit to the vertical profile, although a closer match is 
still desirable. In the channel maps, the data show an 
increasing minor axis thickness with radius, especially in 
the unwarped component, which is not present in this 
model. This suggests a flare. 

Adding a substantial flare to the model (“WF” in fig¬ 
ures) by increasing the scale height with radius allows 
for a better fit to the vertical profile (parameters shown 
in Figure [3]). Additionally, the scale height of the ap¬ 
proaching half decreases slightly for radii larger than 3.1’ 
(13.2 kpc). The emission on the systemic sides of bv di¬ 
agrams is now spread out more, and the aforementioned 
thickening in the channel maps is well reproduced. The 
improvements due to the flare are best seen on the sys¬ 
temic sides (going from systemic to ~75 km s” 1 above 
systemic) of bv diagrams corresponding to radii greater 
than 1.0’ (4.4 kpc) as seen in Figure [5] Note the in¬ 
creased spread between the contours, which is in better 
agreement with the data. Also note the increase in the 
spread along the minor axis. One issue to note is that to 
fit this model, the central scale height is reduced to 3” 
(210 pc), pushing the limits of our ability to constrain 
this parameter given our resolution. However, we do see 
an extremely slight distinction between a scale height of 
3” and 4”, mainly near the peak, shown in Figure [4] In 
the end, the significance of this difference is debatable, 
but we use the 3” (210 pc) scale height for subsequent 
models. 

For the sake of completeness, a second, thicker disk is 
briefly considered, but as one may have already assumed, 
this yields no improvement to the models and is quickly 
dismissed. For instance, the thickening in the minor axis 
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FlG. 2.— Channel maps of NGC 4013 show the substantial warp. Velocities in km s 1 are given in each panel. 
mJy beam -1 , or 6.4*-10 19 cm~ 2 ) and increase by factors of two. The cube is rotatied 66°. 
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FlG. 3.— The parameters used in the final models of NGC 4013 
for the approaching and receding halves. Parameters presented here 
are consistent for all models except for the scale height and inclina¬ 
tion distribution in the case of the 1C model, and the scale height 
distribution for the W model, for which these parameters are held 
constant by definition. Note the relatively symmetric distributions 
of each quantity. The bold line in the bottom-right panel represents 
the inclination 

direction evident in the channel maps as discussed above 
is naturally well matched by a flare. A thick disk would 
not reproduce such behavior and add unnecessary com¬ 
plexity to the models. 


5. THE ADDITION OF A LAG 
Improvements are seen via the addition of a lag of 
—21 tu km s -1 kpc -1 starting at a radius of 1.4’ (5.8 
kpc). The lag causes the observed velocity of the gas to 
be shifted in the direction of the systemic velocity, with 
this shift increasing with height. This causes the ter- 
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FlG. 4.— Vertical profiles of the data, one-component (1C), line 
of sight warp (W), and warp with a flare (WF), summed over a 
range of ±60” (4-3 kpc) from the center in order to avoid as much 
of the warped regions as possible. Additionally, we present the WF 
model with a central exponential scale height of 4 ” as opposed to 
3” in order to illustrate the subtle difference between the two since 
such differences are pushing the limits of the resolution of the data. 
Note how the 4” model slightly underestimates the peak. Given this 
very slight preference, we use the 3” scale height in our final models, 
although the central scale height could indeed be closer to 4 ” • A color 
figure is available in the online version. 


rninal side of the contours in bv diagrams to shift from 
a relatively shallow U-shape to a steeper V-shape. Al¬ 
though we attempt to alter as few parameters as pos¬ 
sible, a slight adjustment to the rotation curve in the 
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FlG. 5.— bv diagrams showing all of the models presented in this paper. Note the clearly poor fit of the 1C model, especially in the narrow 
systemic and wide terminal sides of each panel, as well as near the center. The W model is an improvement, but still lacks the proper 
distribution near the systemic sides that is achieved through the addition of a flare. Contours are as in Figure [M Note the improved match in 
the curvature on the terminal sides of panels corresponding to ± 2.0 ’ and it 2.3 ’ that is due to the lag in this model. 
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Major Axis Offset (ARCMIN) 


FlG. 6.— Representative channel maps showing the approaching 
(A) and receding (B) halves of NGC 4013. Note the splitting in 
the panel corresponding to 736 km s ~ 1 and how this is achieved by 
adding a warp component along the line of sight. The 1C, W, and 
WF models are shown here as these are the models exhibiting the 
most substantial changes from one another. The WFL and WFLv 
models are omitted from this figure as they are very similar to the 
WF model in channel maps since the lag is only in the flat part of 
the disk, which is visually overwhelmed by the outer, warped region 
here. Contours are as in Figure[2\ 



approaching half (190 km s _1 to 200 km s _1 ) yielded a 
better fit for the lag models. The subtle improvements 
from the addition of such a lag are best seen in panels 
corresponding to the terminal sides of ±2.0’-2.3’ of Fig¬ 
ures [5] and [7] (WFL). (Recall that there is a flaring layer 
in NGC 4013, so there is more Hi present at high z at 
large radii, making it easier to detect a lag.) As for the 
central regions within ^1.3’, due to a lack of spatial res¬ 
olution relative to the thickness of the disk, we cannot 
use the bv diagrams to meaningfully constrain the lag in 
those regions. 

We see further improvement by allowing the lag to 
shallow with radius (WFLv) in Figures [5] and m Note 
how the contours on the terminal side in the panels cor¬ 
responding to ±1.7’ and ±2.0’ have a “V” shape on the 
terminal side in both the data and radially varying lag 
models, while those at ±2.3’ are distinctly more rounded 
in both the data and the radially varying lag model than 
in the constant lag model. This behavior can only be 
reproduced if the lag is allowed to shallow radially. The 
shallowing itself is parametrized in Figure [5] 

For completeness, we show lv diagrams of the 1C and 
final WFLv models in Figure [H 


6 . THE UNCERTAINTIES OF THE MODELS 

The rotation curves on each half nearly match each 
other, and the average uncertainty in individual rings is 
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FlG. 7. — Select bv diagrams of the data, warped model with a 
flare (WF), warped model with a flare and constant lag (WFL) , and 
warped model with a flare and radially shallowing lag (WFLv) for 
the approaching (A) and receding (B) halves. The orange lines trace 
the approximate shape of the terminal sides of the bv diagrams in 
the data. In Figure A, note the overall rounded appearance of the 
WF model compared to the data. Note also how the terminal side 
of the WFL model is too pointed in the panel corresponding to 2.3’, 
while the WFLv is somewhat rounded, in closer agreement with the 
data. In Figure B, the WF model is again too rounded compared 
to the data. Regarding the radially shallowing lag, note the rounded 
appearance of the second and third contours on the terminal side of 
the data at —2.3’. The WFL and WFLv models both fit these curves 
quite well at —2.3’, with the WFL model perhaps being slightly too 
pointed. However, at —2.0’ and —1.7’, it becomes evident that the 
lag in the WFL model is too shallow, especially at larger minor axis 
offsets, closer to 950 km s ~ 1 . Thus, to fit the data at smaller radii, 
the lag would need to become steeper for the WFL model, resulting 
in too much of a V-shape in the —2.3’ panel. Contours are as in 
Figure \M A color figure is available in the online version. 



approximately 5 km s -1 . The runs of the position an¬ 
gle and inclination are both extremely sensitive to small 
changes of a degree or two in a single ring. The narrow 
scale height near the center should be considered some¬ 
what skeptically as it approaches the limits of our ability 
to constrain this parameter given the resolution of the 
data and should be considered an upper limit. Still, a 
slight change can be seen in the vertical profile when the 
scale height near the center is changed from 3” to 4” 
(Figure 0]). It is also clear that there is a flaring layer in 
NGC 4013. 

We consider uncertainties in the lag by creating min¬ 
imum and maximum value models which are judged to 
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FlG. 10.— Select bv diagrams of the data, minimum, optimal, 
and maximum lag models for the approaching (A) and receding (B) 
halves. The various lags are parametrized in Figure 0 Note in par¬ 
ticular the subtle differences in the curvature along the terminal sides 
in each model. The minimum lag model is too rounded, while the 
maximum model has a distinct “V” shape not present in the data. 
Contours are as in Figure m 


Fraction of R 25 

FlG. 8.— The distributions of lag magnitudes for the minimum, 
optimal, and maximum lag models shown in Figure ] 1 0\ Constant lag 
values are kept to as large of radii as possible before it is obvious that 
radially shallowing is necessary. Thus, the lag may shallow over its 
entire range. Note the increased uncertainty in the upper limits due 
to the spatial resolution along the minor axis of the data leading to 
smaller changes introduced by steeper lags. 
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FlG. 9.— lv diagrams showing the data, one component (1C), and 
final (WFLv) models. The most substantial changes are near the 
systemic side and at large minor axis offsets. Contours are as in 
Figure \2\ 


be still consistent with the data in addition to the op¬ 
timal lag (Figure [TUI ■ The various lag distributions are 
parametrized in Figured] Note the high maximum lag al¬ 
lowed, which is due to the thinness of the disk in relation 
to the resolution of the data. This difficulty is because of 
the thinning along the minor axis of the terminal side of 
the bv diagrams when steeper lags are introduced. The 
thickness on this side becomes insensitive to the lag for 
relatively steep lags. For a similar reason, an upper limit 
for the magnitude of the lag is increasingly difficult to 
determine with decreasing radius. Within ~7 kpc, our 
small upper limit on the scale height limits our ability to 
constrain steeper lags that result in even further thinning 
along the terminal side of the minor axis. 

By summing the square of the residuals for each model 
interior to 1.15’ (10.2 kpc), we attempt to quantify the 
improvement due to each feature. We do not consider 
the addition of a line of sight warp component (W) to 
a (1C) model here since it is in reality connected to the 
warp component across the line of sight which is obvious 
upon inspection. Adding a flare (WF) to the warped 
(W) model yields a decrease of 3.7% in the approaching 
half and 2.6% in the receding half. The addition of a 
constant lag results in a total increase with respect to the 
W model of 11.3% and a decrease of 4.9%, respectively, 
while allowing the lag to shallow radially provides an 
decrease of 8.1% in the approaching and a decrease of 
5.3% in the receding compared to the W model. 

The rounding on the terminal side of the bv diagram 
at 2.3’ in the approaching half indicates some degree 
of shallowing, and yet the numerical goodness of fit is 
worse than for the model with a constant lag (for the 
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approaching half only). This inconsistency highlights 
the problem inherent in tilted-ring modeling of galax¬ 
ies, in particular their faint and often locally-disrupted 
extra-planar layers, which greatly hinders automated fit¬ 
ting efforts, rendering a by-eye approach necessary as 
demonstrated thr oughout t he liter ature (e.g. Fraternali 
et al. 2004/2005. lOosterloo et all I2Q07L Zschaechner et 
al. 2011/2012/2015. lKamphuis et al.ll2013l [Gentile et all 
2013). The effects of lags predominantly occur in the 
faintest gas - only a fraction of the total emission. Thus, 
the quantifiable improvements from them are suppressed 
over large regions. Furthermore, if bright material close 
to the midplane is displaced by a lag (as is the case here), 
the detriment from that dominates the final quantization 
of the overall improvement to the models. This issue is 
particularly relevant to NGC 4013 given its thin disk, 
which pushes the limits of the resolution of these data. 
However, resolution is not an issue for the fainter, more 
diffuse emission extending above the midplane in which 
the lag is clearly observed and there are indications of its 
shallowing. 

What these models demonstrate is that the overall 
characteristic V-shape of the terminal sides of the bv dia¬ 
grams, as seen in other galaxies with lagging extra-planar 
sc H,i (e.g. NGC 891, NGC 4565, NGC 4244) can only be 
achieved through the addition of a lag. The presence of 
this V-shape is independent of the rotation curve at the 
midplane. Furthermore, the rounding of this V-shape 
at large radii, also demonstrated in the aforementioned 
galaxies, is achieved through shallowing of the lag. It is 
possible that the rounding which we attribute to shallow¬ 
ing here, could be indicative of a constant lag starting 
higher above the midplane with increasing radius (not 
shown). However, such a scenario, in which the round¬ 
ing of the terminal sides of bv diagrams is attributed 
solely to starting height of the lag, is unlikely, although 
it, or a combination of shallowing and variable starting 
height cannot be ruled out based on these data alone. 

7. DISCUSSION & CONCLUSIONS 
7.1. NGC 4013 

The central H I layer in NGC 4013 is rather thin with 
an upper limit of 4” (280 pc) for the central scale height. 
The layer then flares up to 15” (1 kpc) at radii greater 
than 7 kpc. The rotation curve peaks at 195 km s -1 and 
is approximately flat, indicating a large dynamical mass 
(calculated as 1.1 xlO 12 M 0 using a rotational velocity of 
190 km s -1 at a radius of 13.3 kpc). We see no evidence 
for substantial amounts of extra-planar H I, but see evi¬ 
dence for a radially shallowing lag with a peak magnitude 
of—35^28 km s -1 kpc -1 where the Hi layer flares. The 
lag then shallows with radius, going to zero near R 25 . 
The lag is modeled starting from the midplane. 

Our model s agre e remarkably well with those presented 
in IBottemal ( 1996 b the most notable difference being a 
20-30 km s -1 decrease i n rotat ional velocity in the outer 
warped radii in the IBottemal (1996) model, while our 
rotation curve remains flat. This difference is likely due 
to degeneracies in modeling substantially warped regions. 
Recall that the lag is considered only at radii interior to 
the warp, thus this discrepancy has no bearing on the 
final conclusions involving the lag. 

The lack of extra-planar H 1 is intriguing given the 
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FlG. 11.— Measured radial variations of lags in several nearb y 
edge-on spiral galaxies including NGC 891 lOosterloo et a TjUmFl), 
NGC 4%44 1 'Zschaech ner et al\\201 1 \) , NG C 4565 IZs chaechne r et al\ 
1 20121), NGC 5023 (Kanrphuis et al. 20131) , NGC 4302 and NGC 3044 
IZschaechner et al. 2015). Note the general trend of lags beginning 
to shallow near 0.5R 25 , and reaching their shallowest values near 
R 25 itself. A color figure is available in the online version. 

extra-planar DIG and dust. However, it is possible that 
there is a relatively large escape of ionizing radiation 
into the gaseous halo, resulting in prominent extra-planar 
DIG but not HI. 

As a side note, we see no ev idence for the potential H I 
clouds noted bv IBottemal (1 9951 ) that are slightly beyond 
and near the highest point in the warp o n the approach¬ 
ing half and conclude (as IBottemal (|1995lf suggested) that 
they are likely artifacts in t hose data. 

As is also noted bv IBottemal (1199511 . there are no in¬ 
dications of a bar in Hi, although one appears to exist 
in CO and at optical wavelengths dGarcia-Burillo et all 
119991 b The location of the bar they describe coincides 
with the central region of lower column densities before 
the initial peak in Hi (Figure [3j). 

7.2. The Lag in NGC f013 in the Context of Other 
Galaxies 

The lag in NGC 4013 is the steepest measured for H I 
to date. It also shows radial shallowing similar to several 
existing measured lags (Figure ITO) . For a more exten¬ 
sive presentation invol ving lag properties among a larger 
sample of galaxies see IZschaechner et all (2 0151 1. in par¬ 
ticular Table 5. The addition of NGC 4013 does not 
change any conclusions concerning lags and SF proper¬ 
ties already presented in that work. 

7.3. Implications for Theoretical Scenarios 

As has been previously noted, observed lag magni¬ 
tudes are substantially steeper than those predi cted by 
purely ballistic models (e.g. iCollins et al.l 120021) . The 
lag in NGC 4013 is n o exception. Recent models by 
(Marinacc i et al.1 1201 111 consider a momentum exchange 
between galactic fountain clouds and an initially static, 
but subsequently slowly-rotati ng hot halo that is buil t 
up via accretion. As noted by (|Zschaechner et ahll2015jl . 
the radial shallowing observed in these galaxies indicates 
that if lags are indeed due to such a scenario, then this 
momentum exchange must be occurring predominantly 
within R 25 . 

The shallowing itself is not yet understood, and must 
be considered in future simulations. Benjamin (2002, 
2012 ) suggests that extra-planar pressure gradients could 
impact lags and their radial variations. There have 
been few observational constraints to date for such pres¬ 
sure gradients, a situation soon to be remed ied by deep 
continuum observations of edge-on galaxies (llrwin et al .1 
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I2012T ). Their sample includes seven galaxies with mea¬ 
sured lags, including five (as well as the receding side of 
NGC 4302) with radially shallowing lags. 

8 . ACKNOWLEDGMENTS 

We thank the operators at VLA for overseeing the ob¬ 
servations. The National Radio Astronomy Observatory 
is a facility of the National Science Foundation operated 
under cooperative agreement by Associated Universities, 
Inc. This research has made use of the NASA/IPAC 


Infrared Science Archive, which is operated by the Jet 
Propulsion Laboratory, California Institute of Technol¬ 
ogy, under contract with the National Aeronautics and 
Space Administration. We also thank Rene A. M. Wal- 
terbos, Gregory B. Taylor and Peter Zimmer for con¬ 
structive comments. This material is based on work par¬ 
tially supported by the National Science Foundation un¬ 
der grant AST-0908106 to R.J.R. We thank the anony¬ 
mous referee for exceptionally insightful and constructive 
comments regarding the presentation of these models. 


REFERENCES 


Barnabe, M., Ciotti, L., Fraternali, F., Sz Sancisi, R. 2006, A&A, 
446, 61 

Boomsma, R., Oosterloo, T. A., Fraternali, F., van der Hulst, 

J. M., & Sancisi, R. 2005, A&A, 431, 65 
Bottema, R. 1995, A&A, 295, 605 
—. 1996, A&A, 306, 345 

Bottema, R., Shostak, G. S., Sz van der Kruit, P. C. 1987, Nature, 
328, 401 

Collins, J. A., Benjamin, R. A., Sz Rand, R. J. 2002, ApJ, 578, 98 
de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Jr., H. G., Buta, 
R. J., Paturel, G., Sz Fouque, P. 1991, Third Reference 
Catalogue of Bright Galaxies, ed. de Vaucouleurs, G., de 
Vaucouleurs, A., Corwin, H. G., Jr., Buta, R. J., Paturel, G., Sz 
Fouque, P. 

Fraternali, F., Sz Binney, J. J. 2006, MNRAS, 366, 449 
—. 2008, MNRAS, 386, 935 

Garcia-Burillo, S., Combes, F., & Neri, R. 1999, A&A, 343, 740 
Gentile, G., et al. 2013, A&A, 554, A125 

Giuricin, G., Marinoni, C., Ceriani, L., Sz Pisani, A. 2000, ApJ, 
543, 178 

Heald, G., et al. 2011, A&A, 526, A118+ 

Heald, G. H., Rand, R. J., Benjamin, R. A., Sz Bershady, M. A. 
2007, ApJ, 663, 933 

Hodges-Kluck, E. J., Sz Bregman, J. N. 2013, ApJ, 762, 12 
Howk, J. C., Sz Savage, B. D. 1999, AJ, 117, 2077 
Irwin, J., et al. 2012, AJ, 144, 43 
Irwin, J. A. 1994, ApJ, 429, 618 

Irwin, J. A., English, J., Sz Sorathia, B. 1999, AJ, 117, 2102 
Jozsa, G. I. G., Kenn, F., Klein, U., Sz Oosterloo, T. A. 2007, 
A&A, 468, 731 

Kamphuis, P., Peletier, R. F., van der Kruit, P. C., Sz Heald, 

G. H. 2011, MNRAS, 414, 3444 
Kamphuis, P., et al. 2013, MNRAS, 434, 2069 
Leitner, S. N., Sz Kravtsov, A. V. 2011, ApJ, 734, 48 


Li, J.-T., Sz Wang, Q. D. 2013, MNRAS, 435, 3071 
Marinacci, F., Fraternali, F., Nipoti, C., Binney, J., Ciotti, L., & 
Londrillo, P. 2011, MNRAS, 415, 1534 
Martfnez-Delgado, D., Pohlen, M., Gabany, R. J., Majewski, 

S. R., Penarrubia, J., Sz Palma, C. 2009, ApJ, 692, 955 
Miller, S. T., Sz Veilleux, S. 2003, ApJ, 592, 79 
Oosterloo, T., Fraternali, F., Sz Sancisi, R. 2007, AJ, 134, 1019 
Putman, M. E., Peek, J. E. G., Sz Joung, M. R. 2012, ARA&A, 
50, 491 

Rand, R. J. 1996, ApJ, 462, 712 

Rossa, J., & Dettmar, R. 2003, A&A, 406, 493 

Rueff, K. M., Howk, J. C., Pitterle, M., Hirschauer, A. S., Fox, 

A. J., Sz Savage, B. D. 2013, AJ, 145, 62 

Rueff, K. M., Howk, J. C., Wotta, C., Croxall, K. V., Savage, 

B. D., Sz O’Meara, J. 2014, in American Astronomical Society 
Meeting Abstracts, Vol. 223, American Astronomical Society 
Meeting Abstracts #223, #246.19 

Sancisi, R., Fraternali, F., Oosterloo, T., Sz van der Hulst, T. 
2008, A&A Rev., 15, 189 

Springob, C. M., Haynes, M. P., Giovanelli, R., Sz Kent, B. R. 
2005, ApJS, 160, 149 

Tullmann, R., Pietsch, W., Rossa, J., Breitschwerdt, D., Sz 
Dettmar, R. 2006, A&A, 448, 43 

Tully, R. B., Fisher, J. R., Sz Madore, B. F. 1988, JRASC, 82, 305 
Verstappen, J., et al. 2013, A&A, 556, A54 
Wakker, B. P., Sz van Woerden, H. 1997, ARA&A, 35, 217 
Zschaechner, L. K., Rand, R. J., Heald, G. H., Gentile, G., Sz 
Jozsa, G. 2012, ApJ, 760, 37 

Zschaechner, L. K., Rand, R. J., Heald, G. H., Gentile, G., Sz 
Kamphuis, P. 2011, ApJ, 740, 35 
Zschaechner, L. K., Rand, R. J., Sz Walterbos, R. 2015, ApJ, 799, 
61 




